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Underwater  acoustic  systems  detect  the  presence  of  objects  in  the  water  either  by 
directly  sensing  acoustic  energy  radiated  by  the  object  -  passive  sonar,  or  by  transmitting 
an  acoustic  signal  and  detecting  the  reflection  from  the  object  -  active  sonar.  In  general, 
the  acoustic  energy  emitted  by  the  source  reaches  the  receiver  through  different  paths:  as 
direct,  reflected  and/or  diffracted  waves.  The  signals  may  be  received  at  a  single  sensing 
point  or  at  multiple  sensing  points.  The  processing  techniques  depend  on  the  number  of 
receiving  points  and  on  the  number  of  arrivals  at  each  point. 

The  passive  sonar  system  bases  its  detection  and  estimation  of  target  parameters 
on  sounds  which  emanate  from  the  target  itself,  including  machinery  noise,  flow  noise, 
transmission  from  its  active  sonar,  etc.  The  received  signal  represents  the  source  signature 
contaminated  by  noise.  In  addition  to  the  direct  path,  there  are  multipaths  in  sound  wave 
propagation.  The  source  signature  may  include  wideband,  narrowband  and  transient 
components.  These  signals  occur  in  various  multiplicities  and  have  differing  spectra, 
non-planar  wavefronts  and  non-stationary  behavior.  Noise  may  be  partially  correlated 
with  the  signal,  and  high  noise  coherence  may  exist  between  different  sensors. 

In  the  active  sonar  case,  received  echo  represents  the  scattering  signature  of  the  target 
contaminated  by  the  noise.  In  addition  to  self-noise  and  ambient  noise,  the  reverberation 
noise  is  present  and  limits  the  system  performance.  The  reverberation  noise  results  from 
the  scattering  of  the  propagating  signal  energy  due  to  inhomogeneitief  in  the  ocean  and  its 
boundaries.  This  noise  component  is  correlated  with  the  probing  signal  and  proportional 
to  its  power.  Overall,  due  to  relative  target/receiver  motion,  time-varying  signal/noise 
characteristics,  channel  distortion  and  dispersion,  etc.,  the  underwater  signal  detection 
and  estimation  deals  with  non-stationary  signals  and  noise. 

We  were  interested  in  the  class  of  signals  known  as  transients,  which  can  be  described 
as  signals  with  duration  that  is  short  compared  to  the  observation  interval.  They  can  be 
either  deterministic  finite  energy  signals  of  unknown  form,  or  random  finite  duration 
non-stationary  signals.  The  problems  of  detection  and  classification  of  finite  duration 
random  signals,  such  as  the  underwater  acoustic  transients,  have  standard  solutions  based 
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on  the  classical  optimum  detectors  of  random  signals  in  noise[l,  2,  3].  These  detectors 
are  not  designed  to  particularly  distinguish  non-stationary  signals  from  stationary  noise 
background.  However,  it  is  possible  to  take  advantage  of  this  distinguishing  signal  feature 
by  reformulating  the  optimum  detector  in  the  time-frequency  plane.  This  formulation  is 
based  on  the  Wigner-Ville  signal  representation  and  relies  on  the  generalized  singular 
value  decomposition  to  optimally  separate  the  signal  and  noise  subspaces.  This  can  be 
further  applied  to  the  optimal  linear-quadratic  detector  for  problems  involving  certain 
classes  of  non-Gaussian  background  noise[4]. 

Another,  more  fundamental,  aspect  of  this  research  was  the  study  of  probabilistic 
approximation  modelling  of  non-stationary  random  signals.  Recently,  a  new  class  of 
non-stationary  random  signals  has  been  defined,  called  Markov  meshes,  that  is  more 
general  than  the  well-known  Markov  chains  but  with  many  similar  properties  In 
addition,  new  theoretical  results  have  been  obtained  that  enable  us  to  approximately 
model  the  statistical  description  of  a  random  signal  by  simpler,  more  tractable  ones 
(such  as  Markovian  descriptions).  An  objective  of  this  research  was  to  address  the 
approximation  of  random  signals  that  are  given  not  by  their  probability  distributions,  but 
by  their  correlations.  Harmonizable  processes  were  of  particular  interest  here  because  they 
have  well-defined  joint  time-frequency  representations  and  singular  value  spectra.  Various 
classes  of  simple,  tractable  non-stationary  signals  were  analyzed  in  this  framework.  All 
this  allows  for  application  of  the  newly  developed  approximation  theory  to  the  joint  time- 
frequency  representations  and  the  singular  value  spectral  description  of  the  non-stationary 
underwater  acoustic  signals.  This  is  useful  from  several  theoretical  and  practical  points 
of  v  'w.  As  a  result,  approximations  of  joint  representations  and  singular  value  spectra 
of  .  „al  signals  with  the  help  of  the  simpler  ones  are  possible,  providing  at  the  same  time  a 
precise  quantification  of  the  approximation  error.  Then,  the  performance  of  any  detection 
and  classification  scheme  devised  for  the  simpler  types  of  non-stationary  signals  can  be 
easily  evaluated  in  the  presence  of  different  signal  statistics. 
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Signal  Non-Stationarity  as  a  Discriminating  Fepture  for  Detection 
and  Classification  in  the  Time-Frequency  Plane 

We  consider  the  problem  of  detecting  a  zero-mean,  non-stationary,  Gaussian  random 
signal  in  the  zero-mean,  stationary,  colored,  Gaussian  noise  background,  i.e.: 

H0  :  r  =  n 

r,  ,  CD 

ti  \  .  r  —  ti  -f-  .$ 

where  r,  n.  s  are  N-dimensional  vectors  of  samples  representing  the  received  waveform, 
random  noise,  and  random  signal,  respectively.  The  noise  is  assumed  to  have  a  full- 
rank  N  x  :V  covariance  matrix  Rn ,  while  the  signal  A'  x  .X  covariance  matrix  Rs  is 
possibly  rank-defficient.  It  is  further  assumed  that  the  signal  and  noise  are  uncorrelated 
ar.  ,  therefore,  statistically  independent.  Then,  the  detecnon  problem  can  be  rephrased 
in  terms  of  the  following  two  multivariate  normal  models: 

Ho:  r  :  A  (0.  Rn ) 

Hx:  r  :  A  ( 0.  Rn  +  Rs) 

In  words,  the  decision  has  to  be  made  whether  the  received  data  vector  is  a  realization 
of  a  multivariate  zero-mean  normal  distributed  random  vector  with  the  covariance  matrix 
Rn  or  Rn  +  Rs. 

This  problem  has  a  well-known  classical  solution  [1,3]  that  will  be  outlined  next.  We 
shall  then  reformulate  this  solution  in  the  time-frequency  domain,  using  the  Wigner-Ville 
representation,  following  the  approach  in  [5], 

A.  Classical  Formulation 

Standard  solution  to  the  detection  problem  above,  based  on  the  likelihood  ratio, 
compares  the  following  sufficient  test  statistic[l,  3]: 

H  i 

// =  r7^/?"1  -  (Rn  4-  R,)~ 1  r  >  t.  (3) 

Ho 

with  a  threshold  t  determined  according  to  a  suitable  optimization  criterion,  such  as  the 
Neyman-Pearson  criterion,  for  example.  It  is  preferable  to  ’diagonalize’  the  computation 
of  this  statistic  by  means  of  the  generalized  eigendecomposition  of  covariance  matrices 
/?„,  /?,.  Namely,  since  Rn  has  full  rank,  there  exist  a  non-singular  A  x  A'  matrix  V 
such  that[6] 

VrR,V  =  A  =  diag{\,) 

VT  RnV  =  I  =  (hag{\}  (4) 

r  =  ...  r'-v>' 
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where  A,’s  are  the  generalized  eigenvalues  of  R3  with  respect  to  R„,  aranged  in  a  non¬ 
increasing  order  X\  >  Ai  >  . . .  >  Xq  >  0,  Q  =  rank{Rs)  <  X,  and  the  columns  of 
V  are  the  corresponding  generalized  eigenvectors 

Now,  a  diagonalizing  transformation  can  be  applied  to  the  received  data  vector  r 

x  =  VTr  =  [j-i  x>  x3  ...  xy]T  (5) 

and  the  detection  problem  is  transformed  into  a  decision  which  of  the  following  hypothesis 
is  true 

H0:  .r:.V(0.I) 

Hi:  x  :  ;V(0, 1  +  A) 


The  set  of  generalized  eigenvalues,  { A, },  represent  the  distribution  of  signal-to-noise  ratio 
over  the  coordinates  defined  by  the  generalized  eigenvectors.  The  test  statistic  becomes 

Hi 

r/  =  xr  [l  —  (I  +  A)-1]  x  =  xT  [.\(I  +  A ) — 1 1  x  =  xTVx  >  t  (7) 

Ho 

The  resulting  detector  performs  the  linear  transformation,  Eq.  (5),  of  the  received  data, 
with  the  result  being  used  in  the  quadratic  form  computation,  Eq.  (7).  More  explicitly, 
in  terms  of  the  components  of  vector  x  and  the  diagonal  elements  of  f.  A,  the  test 
statistic  is 


Q  Q 

n  =  Z 


A, 

1  +  A, 


Q 


1  +  A, 


According  to  Eq.  (8),  the  optimum  receiver  computes  the  N  squared  correlations  between 
the  received  data  and  the  generalized  eigenvectors,  followed  by  their  weighted  summation. 
In  this  sum,  the  signal  components  along  the  coordinates  with  high  signal-to-noise  ratio 
have  approximately  the  same  weight  l,  while  the  components  along  the  coordinates 
with  low  signal-to-noise  ratio  have  reduced  weight  <  1.  The  test  statistic  >/  is  a  quadratic 
form  in  Gaussian  variates  with  the  conditional  means  and  variances 


E{n  |  Ho)  =  trace) D  =  ~ - 

i=i  1=1  ' ' 

Q 

E(q  |  H\)  —  tr«ce[r(I  +  A)]  =  ^  A, 

i=l 

Q  Q  /  \  V 

Var(r,  \  Ho)  =  2  trace(T2)  =2^  ^  (  fTTV  ) 

Q 

Var{  rj  |  H\)  —  2  trace  T“(  I  +  A)“  =2  ^  A" 

i=l 


(9) 
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The  probabilities  of  detection  and  false  alarm,  Pq  and  Pf\\,  are  governed  by  the 
probability  distribution  of  rj  and  the  threshold  value.  The  distribution  of  the  quadratic 
form  t /  does  not  have  a  neat  analytical  form[2]  and  neither  do  the  probabilities  Pq,  PfA. 
In  such  a  case,  however,  it  is  common  to  use  the  deflection  signal-to-noise  ratio  as  a 
helpfull  indicator  of  the  detector  performance: 


[E(l  |  H\)  —  E(l  ]  //p)]~ 

\'ar(l  |  Ho) 


Q 

2£ 


1=1 


A? 

( 1  +  A,  )* 


(10) 


In  the  high  SNR  case,  most  of  the  generalized  eigenvalues  are  large,  A,  »  1,  and 


D 


Mean  (S.XR)  1 


(ID 


When  the  SNR  is  very  low,  a  situation  known  as  a  low  energy  coherence  case  or  threshold 
detection,  all  generalized  eigenvalues  are  small.  A,  c  1,  and 

Q 

=  ?  Mean  (SNR)2  (12) 

-  1=1 


In  summary,  the  classical  optimum  receiver  consists  of  a  bank  of  correlators,  followed 
by  a  weighted  summer  of  they  squared  outputs  (Eq.  (8)),  and  the  threshold  comparator. 
Its  performance  is  governed  by  the  actual  probability  distributions  of  the  quadratic  form  rj 
under  the  two  hypothesis.  A  gross  indicator  of  the  receiver  performance  is  the  deflection 
signal-to-noise  ratio,  Eq.  (10). 

B.  Time — Frequency  Formulation 

The  optimum  detection  procedure  described  in  the  previous  section  has  an  equivalent 
formulation  in  the  time- frequency  domain.  This  is  of  interest  since  the  signal  to  be 
detected  is  non-stationary,  a  problem  for  which  intuition  suggests  comparison  of  time- 
frequency  signal  representations  [5]. 

The  time-frequency  domain  equivalents  of  the  inner  products  in  the  classical  optimum 
receiver,  Eq.  (8),  can  be  obtained  using  any  transformation  that  preserves  the  inner  prod¬ 
ucts.  Among  such  transformations,  Wigner-Ville  representation!?]  has  certain  advantages 
that  make  it  the  most  attractive  for  our  purposes.  Namely,  it  is  the  only  one  that  has  all 
of  the  following  desirable  properties[7,  8]:  (1)  it  is  real-valued;  (2)  it  is  invariant  to  nme 
and  frequency  shifts;  (3)  it  has  the  same  region  of  support  as  the  signal;  (4)  it  preserves 
the  inner  products;  (5)  it  has  the  proper  marginals;  (6)  among  all  the  representations 
satisfying  (1) — (5)  it  is  the  one  best  localized  in  the  time-frequency  plane. 

It  should  be  noted,  however,  that  these  properties,  while  valid  in  the  continuous-time 
case,  are  not  generally  valid  for  the  discrete-time  Wigner-Ville  representation  which  is 
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periodic  in  frequency  with  period  k  [7].  All  of  the  above  properties  are  valid  only 
if  the  discrete-time  signal  spectrum  is  non-zero  over  an  interval  less  than  tc  on  its 
fundamental  period  2t  [7],  Two  practically  important  such  cases  are:  (a)  when  the 
real-valued  continous-time  signal  is  oversampled  by  at  least  a  factor  of  2,  and  (b)  when 
the  signal  is  analytic,  i.e.  its  spectrum  vanishes  for  the  ’negative’  frequency  half  of  the 
fundamental  period  2 it  [7].  In  the  rest  of  this  papier,  it  will  be  assumed  that  the  signals 
of  interest  satisfy  the  aforementioned  condition. 

a.  Wigner-Ville  Representation  of  Discrete-Time  Signals 

Let  x{k),  y{k).  x.  y  6  C.  k  €  Z  be  complex-valued  discrete — time  signals,  for 
which  there  exist  Fourier  transforms  X(0).  V(0).  Then,  the  discrete-time  cross- Wigner- 
Ville  representation  of  x  and  y  is  defined  by[7] 

\VIy(k.0)  =  2  Y.  -rU* +  03) 

m  =  —  y c 

Here,  can  be  interpreted  as  the  Founer  transform  of  the  sequence  cfc(m)  = 

x(k  +  m)  y*(k  -  w)  considered  as  a  function  of  m  for  fixed  k.  Obviously,  the  cross- 
Wigner-Ville  representation  is  complex-valued  and  penodic  in  frequency  with  penod 
t,  and  \Viy(k.O)  =  \V*r{k.O).  It  can  be  expressed  in  terms  of  the  signal  Founer 
transforms: 

7T 

Wiy(k\0)  =  ~  f  X(0  +  Z)Y'i0-‘)(>*k  (14) 

—  7T 

Of  particular  interest,  however,  is  the  special  case  y  =  x  that  gives  the  plain  discrete¬ 
time  Wigner-Ville  representation: 


WAk.B)  =  2  Y  x(A’  ^  v"'  -r*ik  ~ 

ITl  =  —  X. 

T 

=  i  J  X(0  +  £)  X*if)-Z)<J*k  <lz 


(15) 


which  is  real-valued.  I Vr{k,0)  =  \V*(k.O),  and  symmetric  with  respect  to  frequency, 
Wx(k.B)  -  l Vr*{k.  -0).  For  the  detection  problems,  the  most  important  property  of  the 
Wigner-Ville  representations  is  expressed  by  the  Moyal’s  formula,  which  in  discrete-time 
case  reads[7] 


xz  f  Y  ii :ririik.0)\v;^L0uiH 

-in 

=  (xj.  X  f)(x2.  X4)‘  +  ( .r  1 .  MyX  .  .Vl-.r.,  |* 

=  (A  i .  A  i )( A  A  i )  4-  ( A  i .  <ST  A  t ;( A  >T  A  1 1 


(16) 
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where 


•*4  1  f 

(xi , x> )  =  V]  Xi {k)  x*(k),  ( A" i ,  A 2 )  =  ~ —  Xi{0)  X2{0)  (W 

- - - -  “'4  (1?) 

.Vf ;rX  =  ejkirx{k)  =  (  — l)*x(A-),  5tA  =  A(0  -  tr) 

However,  assuming  that  the  signals  of  interest  have  the  frequency  domain  support  less 
than  7r,  the  regions  of  support  of  X\(9)  and  Sv  Aj  =  A3 ( 0  -  ?r)  do  not  overlap,  aliasing 
is  avoided,  and  Moyal’s  formula  reduces  to2 


l 

2x 


\VXll2(k.0)\V'IiTi(k.0)d0 


-  Ui--r.j)(.r_>.x4r  =  (-Vi,  .Ai)(Aj.  A.(  )* 


(18) 


For  example,  this  form  of  the  Moyal’s  formula  is  applicable  to  real-valued  signals 
oversampled  by  at  least  a  factor  of  2,  and  to  analytic  signals,  where  the  limits  of 
integration  should  be  taken  to  be  0,  x 

If  the  signals  of  interest  have  finite  discrete-time  support  region  T  of  length  A  = 
2/\  +  1,  i.e.  if 


x(k)  =  y(k)  =  0.  kg  T 


(19) 


then  it  follows  that 


\Vty(kJ))  =  0,  kg  T 


(20) 


and 


Cfc(m)  =  x(k  +  m)  ij*(k  —  rn)  =  0.  |m|  >  I\ 


(21) 


Since  Wty{k.  |)  is  a  Founer  transform  of  the  sequence  cjt(m)  =  x(k  +  rn)  <j*(k  -  rn) 
that  has  finite  duration  m  this  case,  for  each  k,  WT!l(k.O)  is  completely  specified  by 
V/  >  V  frequency  samples  taken  over  its  fundamental  period  rr.  These  samples  are 
given  by  the  discrete  cross-Wigner-Ville  representation 


\Vry(kJ)  =  Wry(k.0=^jl) 

K 

=  2  x{k  +  m)  y*[k  —  m)  . 

m  =  —  h 

=  2  VFT{zk{rn)} 


(22) 


where  k,  /  are  the  time-index  and  the  frequency-index,  and  T.  Q  are  the  discrete¬ 
time  and  discrete-frequency  support  regions,  respectively.  For  the  real-valued  signals 


-  In  Ref  [9]  it  is  incorrectly  stated  that  oversampling  can  not  eliminate  the  aliasing  terms  in  the 
Moyal’s  formula. 


{}  =  {/  |  -  L  <  l  <  L,  M  =  2L  +  1},  and  for  the  analytic  signals  il  =  {/  j  ()</<.)/}. 
Setting  y  —  x  gives  the  plain,  real-valued,  discrete  Wigner-Ville  representation 


K 

Wt(kJ)  =  2  Y  x[k  +  w)x*(k-m)f-^lw. 

rn  =  —  l\ 


k  €  T 


/  €  <> 


(23) 


The  discrete-time  discrete-frequency  version  of  the  Moyal’s  formula  is  obtained,  in 
the  absence  of  aliasing,  by  discretizanon  of  Eq.  (18), 


‘  tell  k£ T 

=  (•*■!. a-;i)(x>,.r4)"  =  (A'i,  A.i)(.Y.>.  .V4)* 


(24) 


This  is  the  equation  that  allows  easy  reformulation  of  the  classical  optimum  detection 
procedure  in  the  rime-frequency  plane.  However,  it  is  necessary  to  consider  the  Wigner- 
Ville  representation  of  random  signals[  10]  first. 

b.  Wigner-Ville  Spectrum  of  Random  Signals 

Let  x.  y  be  the  realizations  of  harmonizable  zero-mean  complex-valued  discrete-time 
random  processes  .V,  Then,  the  following  expectation  exists[10] 

n'xy(k-.e)  =  E{\vxy(k-.o)} 

DC 

=  2  V  E{x(k  +  /»);/*( k  -  m)} 

rr,  =  -x  (25) 

=  2  Y  l?r,(k-  +  >n.k--m)<-J-t)r" 


and  is  called  the  cross-Wigner-VUle  spectrum  of  processes  .V.  As  its  deterministic 
counterpart,  it  is  periodic  with  period  We  see  that  VV  (k.  - )  is  the  Founer 
transformation  of  the  cross-covanance  RTy(k  +  m.  k  -  m )  with  respect  to  displacement 
m.  Setting  y  =  -V,  we  get  the  plain  real-valued  Wigner-Ville  spectrum  of  ,V  [10] 

DC 

n'x(kj))  =  E{\Vr(k.0)}  =  2  Y  R^k  +  -  '")  '  ~J (26) 

rn  =  —  dc 


which  has  the  same  periodicity,  symmetry  and  support  properties  as  its  deterministic 
counterpart.  In  particular,  if  the  random  signal  is  analytic  (having  the  power  spectrum  that 
vanishes  on  the  ’negative’  frequency  half  of  the  fundamental  period  2rr),  or  oversampled 
by  at  least  a  factor  of  2.  aliasing  along  the  frequency  axis  is  zero. 

If  the  random  signal  is  stationary,  Rr(k  +  w.  k  -  w)  =  Rr(2m  i.  and  the  Wigner- 
Ville  spectrum  reduces  to  the  aliased  ordinary  power  density  spectrum 

X 

W. v ( k.  6)  =  2  Y  Rr(2w)(  ->im*  =  {())  +  Sr  ( 0-  ~  \  (27) 

m  =  — dc 
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However,  if  the  random  signal  is  adequately  sampled,  Sz{9)  and  Sz(0  —  -)  do  not 
overlap,  and 

Wx(k.B)  =  SAB),  OZ®  (28) 

where  $  is  an  interval  of  length  7 r  that  includes  the  region  of  spec-ral  support  as  its 
proper  subset. 

If  the  random  process  .V  has  finite  discrete-rime  region  of  support  T  of  length 
.V  =  2  A'  +  1,  then 


W'x(h.O)  =  0.  k$  T 


(29) 


and 


Rz(k  +  m ,  k  —  in)  —  0,  \m\  >  I\ 


(30) 


For  each  k,  then,  W,y(A\  9)  is  completely  specified  by  M  >  V  frequency  domain  samples 
taken  over  its  fundamental  period  tt.  The  result  is  known  as  the  discrete  Wigner-Ville 
spectrum 

Wx(kJ)  =  Wx(k,9  =  -^/) 

K  p  ^  J 

=  2  £  Rz(kx  m.k  -  in) 
m  =  —  A’ 


(31) 


/  €  n 


=  £’{Hx(A-./)} 

where  T.  H  are  the  discrete-time  and  discrete-frequency  support  regions,  respectively. 

Infinite  duration  random  processes  are  usually  observed  over  finite  intervals,  and  their 
covariances  may  be  known  or  estimated  only  on  a  grid  of  limited  extent.  In  such  a  case, 
it  is  possible  to  obtain  an  approximation  of  the  sampled  Wigner-Ville  spectrum,  smeared 
one-dimensionally  in  the  direction  of  frequency,  by  means  of  the  discrete  Wigner-Ville 
spectrum,  Eq.  (31). 

For  example,  let  .V  be  a  harmonizable  infinite  duration  discrete-time  random  process 
observed  on  the  interval  —IK  <  k  <  2I\ ,  with  spectral  support  less  then  t  and  the 
covariance  RT(p,q)  known  on  the  gnd  —  2/V  <  p.q  <  2K.  Then.  Rz{k  4-  m.k  -  m)  is 
completely  known  on  the  gnd  -K  <  k.w  <  K  and  the  sampled  approximation  of  the 
Wigner-Ville  spectrum  can  be  computed  as 


2C 

MTv(t,6l  =  -^/)  =  2  Y,  g(m)Rr(k  +  m.k-m)c-J$lm 

m=  —  oc 

=  \G(0)  *0  Wx{k.9)\n=fji 
K 

=  2  Rx(k  +  rn,k  —  nt) 

in  =  —  A 

=  Wx(k.l).  I  e  <>.  - K  <  k  <  K 


(32) 


where  *0  denotes  convolution  along  the  9  axis,  </(  m )  is  the  rectangular  truncating  window 
and  G'(^)  is  its  Founer  transform.  If  the  process  .V  is  also  stationary,  since  there  is  no 
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aliasing  in  the  time-frequency  plane,  the  Wigner-Ville  spectrum  equals  the  ordinary  power 
spectrum,  Eq.  (28),  and  its  sampled  smeared  approximation  is 

Wx(kJ)  =  [G{9)  *  S_rlt))}d=±l  =  £(/),  lea  1A-I  <  A'  (33) 


In  words,  the  finite-record  approximation  of  the  Wigner-Ville  spectrum  of  an  adequately 
oversampled  stationary  process  is  equal  to  the  finite-record  approximation  of  its  ordinary 
power  spectral  density. 


c.  Detection  Based  on  the  Discrete  Wigner-Ville  Representation 

The  detection  problem,  Eq.  (1),  can  be  rephrased  in  the  time-frequency  plane  as  the 
choice  between  the  hypotheses 


H0  :  U  V  =  Wn 
H{  :  ll'r  =  II  „  +  \VS 


(34) 


where  1V„.  U-3  are  the  discrete  Wigner-Ville  representations  of  the  particular  noise  and 
signal  realizations. 

We  are  ready  now  to  transform  the  time-domain  correlator  based  optimum  receiver, 
Eq.  (8),  into  the  time-frequency  domain  correlator  based  receiver.  The  inner  products  in 
Eq.  (8)  can  be  replaced,  according  to  the  discrete  Moyal’s  formula,  Eq.  (24),  by 


r 


r 


1 

lu 


YYU 7(M  w\(k.l) 

t€G  f:6 T 


(35) 


where  \Yr{k.  /)  is  the  Wigner-Ville  representation  of  the  received  data  vector  and  W,(k\  I) 
are  the  Wigner-Ville  representations  of  the  generalized  eigenvectors  v{l) .  Using  this  result 
in  Eq.  (8),  we  get 


// 1 

>  t 

leu  ke T 

no 


where  B(k,l)  is  the  weighted  sum 


Q 


Ba-l]  =  :,uY. 


A, 


23/  4— »  l  +  A, 

1=  1 


U \(k.l) 


(36) 


(37) 


which  is  fixed  by  the  generalized  eigendecomposition.  Eq.  (4),  does  not  depend  on 
the  data  vector  r,  and  can  be  precomputed.  If  the  slices  of  U'r  and  B  for  the  fixed 
time-index  k  are  labeled  a-J  =  [I Vr ( A. 0 )  Wr{k,  1)  ...  Wr\k..\I  -  l)]  and  h[  = 
[B(k.  0)  B(k.  1)  ...  B(k.  M  -  1 )],  respectively,  the  optimum  test  has  a  more  explicit 
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Y,  hk  '  U'* 

keT 


> 

< 

Hi 


(38) 
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Accordingly,  the  time-frequency  domain  equivalent  of  the  classical  optimum  receiver 
computes  the  slice  by  slice  cross-correlation  of  the  received  data  Wigner-Ville  repre¬ 
sentation  with  the  precomputed  weighted  sum  of  Wigner-Ville  representations  of  the 
covanance  matrix  generalized  eigenvectors.  An  alternative,  expression  for  the  optimum 
test  is  obtained  by  introducing  the  matrix  notation  !l'r  =  [irr(A\/)j.  H  =  {B{k.h\ 

H  i 

ij  =  trace  j  B  U"/  |  =  trace  ^BT  llr|  >  t  (39) 


//,) 

In  the  reference  [5],  a  similar  reformulation  of  the  continuous-time  optimum  detector  has 
been  presented  for  the  case  of  white  noise  only. 

The  test  statisnc  conditional  expectations  remain  the  same  as  those  of  the  classical 
receiver,  as  given  by  Eos.  (9).  It  is  interesting  to  relate  those  results  with  the  Wigner- 
Ville  spectra  of  signal  and  noise.  Taking  the  conditional  expectations  of  Eq.  (39),  using 
the  Eqs.  (31)  and  (32),  and  combining  the  results  with  Eqs.  (9)  we  get 

Q 


E('1  I  Ho)  =  tracers  VV'[-j  = 

i=i 


A, 


1  -r  A, 


E(<] 

1  //.  )  = 

trace 

i*( 

vv[-  +  vvj 

)  -E 

_ _ 1 

■  .  — _  .  >1 

i —  I 

Q 

Var{ 

>/ !  //o) 

-  '2  trace 

t 

=  -’E 

i=i 

\'ar{ 

n  I  H 1) 

—  2  trace < 

[[«(vv(-  +  wj)]'' 

A. 


(40) 


where  VV’v  is  the  finite-record  approximation  of  the  Wigner-Ville  spectrum  of  noise,  as 
discussed  in  connection  with  Eq.  (32),  and  VV5  is  the  discrete  Wigner-Ville  spectrum  of 
a  finite  duration  non-stationary  signal,  as  in  Eq.  (31).  The  deflection  signal-to-noise  ratio 
is  the  same  as  before  and  given  by  Eqs.  (10)  —  (12). 

We  see  that  an  alternative  time-frequency  domain  realization  of  the  optimum  receiver 
is  possible.  It  matches  optimality  with  physical  interpretation  provided  by  the  use 
of  the  Wigner-Ville  representation,  which  is  particularly  relevant  for  the  detection  of 
non-stationary  signals.  However,  the  receiver  performance  remains  the  same.  The  real 
significance  of  this  time-frequency  formulation  is  that  it  allows  incorporation  of  additional 
information  that  will  bring  improvement  in  performance,  which  is  the  subject  of  this 
proposal. 


C.  Subspace  Based  Detection 

Our  objective  is  to  use  a  prion  available  information  that  signal,  as  opposed  to  the 
noise,  is  non-stationary  in  order  to  suppress  that  background  noise  and  thereby  improve 
the  deflection  signal-to-noise  ratio  and  detector  performance.  How  do  we  accomplish 
this  is  discussed  next. 


1 1 


a.  Optimum  Signal/Noise  Subspace  Decomposition 

The  Wigner-Ville  spectrum  of  the  stationary  random  signal,  Eq.  (28),  and  its  finite 
record  approximation,  Eq.  (33),  are  simple  outer  products  of  the  power  spectral  density 
as  a  function  of  frequency  and  a  constant  function  of  time.  Consequently,  the  matnx 
of  samples  of  the  finite  record  approximation  of  the  Wigner-Ville  spectrum  of  stationary 
noise,  W.v,  is  a  unit-rank  matrix.  On  the  other  hand,  the  discrete  Wigner-Ville  spectrum 
of  the  random  acoustic  transient,  Eq.  (31),  represented  by  the  matrix  W’s,  has  rank  that  is 
generally  much  higher  than  one.  This  distinction  can  be  exploited  to  suppress  stationary 
noise  and  improve  the  detector  performance. 

The  generalized  singular  value  decomposition  of  the  signal  Wigner-Ville  spectrum 
with  respect  to  the  noise  Wigner-Ville  spectrum  provides  a  necessary  tool  to  accomplish 
this.  Assuming  that  VVy  and  W5  are  A  x  37,  37  >  .X  matrices,  then  there  exist  a 

non-singular  .V  x  .V  matrix  A'  and  orthogonal  matrices  73(37  x  37),  V{.\I  x  37)  such 

that!  11] 

W5  -  A'-1  DsUt 
W.v  =  X~lDnVT 

DS(.X  x  37)  =  diag(ai).  0  <  oi  <  qi  <  ...  <  a y  (41) 

D, j(.V  x  37)  =  diag(3l),  3X  >  ...  >  3q  >  3q+x  =  ...  =  3\  =  0. 

q  =  ranAr(VVv') 

where  x(0,  the  rows  of  .V,  are  the  generalized  singular  vectors  of  W\'  and  VV5  satisfying 

3;  W5WJ  xU]  =  a;  ,r(,)  (42) 


Note  that  for  the  stationary  noise  rank(Wx)  =  q  —  1  and,  therefore,  only  3 1  7^  0.  Let’s 
introduce  the  following  matrix  partitions 


.A  = 


rur 

A’. 


,  C  =  [u(1)  i'll  V  =  [C1*  Vi],  D* 


oi  0 

0  73 1 


(43) 


where 

A3  =  A3  (A'  —  1  x  A’),  D\  =  Di{.X  —  1  x  37  —  1)  =  diag(a> . ay) 

(3  =  (3(37  x  37  -  1),  Vi  =  V] (37  x  37  -1)  (44) 


Then, 


XXVS  =  DSL 


A' W.v'  -  DnVl  = 


where  Oi  is  a  .V  — 


.f  |  Ql  u* 1 

DiUf 

0, 

1  x  37  zero  matrix.  More  specifically 


A3 VVs  =  Dxl  [  =  Y 
•V,  WA-  =  0, 


(45) 


(46) 
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i.e.,  A'i — transformation  maps  the  Wigner-Ville  spectrum  of  stationary  colored  noise  into 
an  all-zero  matrix  Oi.  At  the  same  time,  since  <  Q;  is  small,  the  signal  Wigner-Ville 
spectrum  is  mapped  with  the  minimum  loss  of  signal  power. 

We  can  interpret  these  results  as  saying  that  all  of  the  noise  energy,  together  with 
a  minimal  fraction  of  signal  energy,  lie  along  the  first  generalized  singular  vector  .r(1), 
and  there  is  no  noise  energy  in  the  subspace  determined  by  the  remaining  generalized 
singular  vectors  that  constitute  the  rows  of  A'i[ll]. 

b.  Detection  with  Noise  Subspace  Elimination 

The  generalized  singular  value  decomposition  of  VVy  and  W^,  Eq  (41),  thus  gener¬ 
ates  an  optimum  stationary/non-stationary  subspace  decomposition,  Eq  (45),  that  will  be 
the  basis  of  the  new  detection  scheme  discussed  next. 

Rather  than  correlating  the  received  data  Wigner-Ville  representation  with  a  reference, 
as  in  Eq.  (39),  we  propose  to  first  use  the  A'i — subspace  mapping.  It  is  easy  to  show  that 

H0:  E(X[\Vr)  =  A’lWy  =  0, 

Hi  :  E(Xi\Vr)  =  A',Wa-  +  A'j W5  =  Y 

Therefore  we  propose  to  compute  the  test  statistic 

// 1 

7  =  trace j(A'i£)r(A'iirr)j  <  /i  (48) 

Ho 


where  t 1  is  the  threshold  determined  by  the  suitable  criterion,  such  as  Neyman-Pearson’s, 
for  example.  With  notation  B\  =  A'i  B,  the  conditional  means  can  be  expressed  as 


E{rj\Ho)  =  trace  B[ ( X i VV  ) 


0 


E(r]\H\ )  =  trace  Bj ( A' i  VV\'  +  A'lW's) 


trace  ( B 


<lr) 


(49) 


The  conditional  variances  cannot  be  brought  into  a  neat  analytic  form  but,  intuitively, 
we  expect  them  to  be  reduced  with  respect  to  the  classical  case  due  to  the  elimination 
of  noise  subspace  by  the  A'i — subspace  mapping.  This  point  needs  to  be  verified  by 
Monte-Carlo  simulations,  as  will  be  reported  in  a  future  publication. 


Probabilistic  Approximation  Modelling  of  Non-Stationary 
Underwater  Acoustic  Signals 

Given  an  arbitrary  random  signal  £,  let  A  be  some  class  of  random  signals  with  the 
help  of  which  it  is  reasonable  to  try  to  approximate  £  (e.g.  having  the  same  sets  of 
states  as  <f,  a.s.o,  and  some  sufficient  tractable  structure  of  probability  distributions). 
A  method  has  been  developed  by  M.  Rosenblatt-Roth,  a  consultant  on  this  project, 
permitting  the  explicit  determination  of  (1)  the  random  signal  7  -  £1  €  .1  which  is 
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the  best  approximation  of  £  in  .4,  and  (2)  the  least  error  r  i(£)  corresponding  to  this 
approximation.  This  approximation  has  been  made  in  a  meaningful  sense  so  that  this 
closeness  implies  closeness  in  variation.  Taking  for  .4  various  classes  .4l  kK  1  <  k  <  ii. 
let  the  best  approximants  of  £  be  £*,  and  let  the  corresponding  errors  be  rk(^ ).  1  <  k  <  n. 
If  r/( £)  is  the  smallest  of  them,  we  should  model  the  given  signal  (  by  ^  G  ,4(/).  For 
this  modeling  to  be  effective,  it  is  necessary  (1)  to  chose  the  classes  ,4(i)  in  an  adequate 
manner,  (2)  to  take  n  sufficiently  large,  and  (3)  to  verify  that  the  main  properties  for 
which  we  study  the  signal  £  are  not  lost  in  the  modeling. 

Among  the  best  known  and  most  analyzed  non-stationary  random  signals  are  those 
modelled  by  Markov  chains.  They  have  been  studied  since  the  beginning  of  the  century, 
but  it  was  not  until  1950’s  that  a  significant  progress  was  made.  At  that  time,  R.L. 
Dobrushin[l2]  obtained  various  important  results  concerning  the  central  limit  theorem 
for  Markov  chains.  M.  Rosenblatt-Roth  has  investigated  various  aspects  of  the  law  of 
large  numbers  for  sequences  of  random  variables  forming  Markov  chains  in  the  context 
of  a  study  of  the  transmission  of  signals  produced  by  non-stationary  sources  through  non- 
stationary  channels[13,  14,  15,  16].  He  also  obtained  some  results  concerning  problems 
of  stability  and  compression. 

Methods  will  be  devised  for  the  reduction,  in  non-trivial  ways,  of  more  complex 
non-stationary  random  signals  onto  Markovian  ones.  Therein  lies  the  significance  of 
the  non-stationary  Markov  chains  in  this  research.  For  this  reason,  on  the  basis  of  the 
well-known  results,  new  aspects  concerning  the  non-stationary  Markov  chains  will  be 
studied.  The  basic  idea  of  Markov  chains  was  not  exhausted  by  their  strict  definition.  It 
is  possible  to  extend  this  concept  by  interpreting  the  Markov  idea,  from  a  somewhat  more 
general  point  of  view,  that  the  "future"  of  a  Markov  random  process  is  "independent" 
from  its  "past"  if  its  "present"  is  known.  For  Markov  chains,  the  "present"  is  represented 
by  a  given  finite  set  of  random  variables  immediately  preceding  the  future.  By  relaxing 
the  above  condition  so  that  the  "present"  be  represented  by  a  finite,  but  variable  number 
of  random  variables  preceding  the  "future",  M.  Rosenblatt-Roth  has  recently  obtained  a 
much  larger  class  of  random  signals  than  the  Markovian  ones[17].  This  new  class  of  non- 
stationary  random  processes,  called  Markov  meshes,  can  be  of  Gaussian  or  non-Gaussian 
type  and  has  important  similarities  with  Markov  chains. 

More  precisely,  let  £*  (1  <  k  <  n)  be  a  sequence  of  random  variables;  we  say  that  it 
is  a  Markov  chain  of  order  r  >  1  if,  for  any  k  (r  <  k  <  n),  the  "past"  represented 
by  the  random  variables  —  sk-r-\}  is  independent  of  the  random  variable 

.ft  representing  the  "future"  if  it  is  known  that  the  random  variables  {^._r . *k-\} 

representing  the  "present"  take  given  values  {.r*_r . r*_i}. 

On  the  other  hand,  a  one-dimensional  Markov  mesh  is  defined  in  the  following  way 
by  M.  Rosenblatt-Roth.  For  the  same  sequence  of  random  variables  as  above,  consider 
a  subset  Sk  of  the  index  set  I  —  {1,2.  ....k  —  1}  with  Sk  denoting  the  complement  of 
Sk  with  respect  to  /,  and  let  £(Sk)  =  {£,,  j  e  Sk},  t(Si)  =  {sc7,  j  €  Sk}.  We  say 
that  the  sequence  ( 1  <  k  <  n)  is  a  one-dimensional  Markov  random  mesh  if  for  any 
k  (1  <  k  <  n),  the  "past"  represented  by  random  variables  i{Sk)  is  independent  of  the 
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random  variable  represenring  the  "future"  if  it  is  known  that  the  random  variables 
£(5*)  representing  the  "present"  take  given  values  { ij .  j  e  S^}. 

The  class  of  one-dimensional  Markov  meshes  not  only  contains  the  Markov  chains 
as  a  sub-class,  but  shares  with  them  some  important  properties.  For  this  reason.  Markov 
meshes  appear  to  have  promising  future  in  signal  theory.  The  concept  of  a  Markov  chain 
was  extended  to  multi-dimensional  signals  under  the  name  of  Maikov  field.  However, 
because  of  the  topological  differences  between  a  one-dimensional  sequence  of  sites  and  a 
multi-dimensional  array  of  sites,  this  concept  remains  still  far  from  that  of  a  Markov  chain, 
lacking  many  of  its  essential  properties.  For  this  reason,  in  many  important  applications 
Markov  fields  are  unfit  as  a  model  of  a  multi-dimensional  random  signal. 

M.  Rosenblatt-Roth  has  recently  shown  that  the  basic  idea  of  a  Markov  chain  is  not 
exhausted  by  this  extension.  He  has  extended  the  concept  of  a  one-dimensional  Maikov 
random  mesh  to  a  multi-dimensional  mesh  retaining  the  basic  Markovian  idea  that  the 
"future"  of  the  signal  is  independent  of  its  "past"  if  its  "present"  is  known[17].  In  this  way, 
we  obtain  a  very  large  class  of  non-stationary  random  signals,  which  can  be  of  Gaussian  or 
non-Gaussian  nature,  and  which  posses  important  properties  similar  to  those  of  Markov 
chains.  More  than  that,  the  multi-dimensional  Markov  meshes  admit  one-dimensional 
causal  representations,  a  property  that  reduces  the  study  of  multi-dimensional  signals 
to  that  of  corresponding  one-dimensional  ones.  Among  the  multi-dimensional  random 
signals,  the  multi-dimensional  Markov  meshes  play  the  same  role  as  do  the  Maikov 
chains  among  the  one-dimensional  ones. 

Considering  the  complexity  of  non-stationary  non-Gaussian  random  signals,  an  at¬ 
tractive  manner  of  handling  them  would  be  by  approximating  their  statistical  descriptions 
with  the  help  of  various  adequate  classes  of  simpler  ones.  M.  Rosenblatt-Roth  defined  and 
studied  a  manner  of  approximation  that  allows  the  explicit  determination  of;  (a)  the  best 
approximant  in  the  given  class  of  statistical  models,  and  (b)  the  least  error  committed  in 
this  best  approximation[17,  18,  19,  20]. 

Let:  (.V.  5)  be  a  measurable  space,  L  =  L(X.S)  the  convex  set  of  probability 
measures  on  (X.  5);  A  an  idempotent  operator,  non-linear  in  general,  defined  on  L  with 
values  in  IJ  =  A(L)  C  L.  If  l ,  V  £  £.,  let  h(P  \Y)  be  the  relanve  entropy  of  U 
with  respect  to  V.  If  P  €  L  and  AP  €  L'  is  absolutely  continuous  with  respect  to 
Q  €  V  (AP  <  Q),  let 

h(AP  :  Q:  P)  =  j  P(dX)  log[(  AP)(clX  )/Q(dX )]  (50) 

V 

A  is  the  projector  on  V  and  AP  €  V  is  the  .1-projection  of  P  £  L  on  L'\  A  is  a  regular 
projector  if  for  any  P  €  L, 

rnin{h(AP  :  Q\  P)\  Q  G  L' }  =  0  (51) 

with  equality  iff  Q  =  AP. 

The  quantity  I±(P )  =  h(P  :  AP)  is  the  amount  of  A-information  determined  by 
and  the  quantity  I\  =  h(AP  :  P)  is  the  amount  of  conjugate  A-information 
determined  by  P  e  L.  M.  Rosenblatt-Roth  has  obtained  the  following  results: 
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(1)  If  P  e  L,  then  I  A{P)  >  0,  IA(P)  >  0  with  equalities  iff  P  =  AP,  i.e.  iff  P  6  L' . 

(2)  If  A  is  a  regular  projector  and  IA(P)  is  finite,  then 

nun{h(P  :Q)\  Q  £  L1}  =  lA[P)  (52) 


and  this  minimum  is  reached  iff  =  AP. 

(3)  Let  P(t)  be  a  Founer  transform  of  P  €  L  and  j|.||  the  total  variation  of  a 
completely  additive  function  defined  on  (X.S).  Then 


sup 


{!™ 


lAP)(t) 


\P-AP\\~  <2  wm{IA{P).I\(P)} 


(53) 


From  (1)  it  is  seen  that  IA(  P)  is  good  as  a  measure  of  the  deviation  of  P  from  =  .1(  L ); 
from  (2)  that  AP  is  the  best  approximant  of  P  in  A(L)  and  from  (3)  that  if  IA(P)  is 
small,  so  is  P  -  AP  in  total  variation  and  AP  is  uniformly  close  to  P. 

Considenng  vanous  particular  classes  of  one-dimensional  and  multi-dimensional 
signals,  M.  Rosenblatt-Roth  calculated  the  explicit  expressions  for  the  /i-information 
using  the  expression 

lA(P)  =  f  P(dX)  log[P(d.\)/(AP)(dX)}  (54) 


which  can  be  finite  only  if  AP  <  P. 

Over  the  past  three  years,  M.  Rosenblatt-Roth  has  extensively  studied  such  approxi¬ 
mation  of  one-dimensional  random  signals  with  (a)  sequences  of  independent  identically 
distributed  random  variables,  (b)  non-Gaussian  Markov  sequences  of  order  r  >  1,  (c) 
non-stationary  Gaussian  sequences,  (d)  non-stationary  Gauss-Markov  sequences  of  or¬ 
der  r  >  1,  (e)  non-stationary  Gauss-Markov  meshes,  (f)  stationary  Maikov  chains  of 
order  r  >  l,  (g)  stationary  Markov  meshes  of  vanous  classes.  For  multi-dimensional 
random  signals  M.  Rosenblatt-Roth  has  studied  the  approximation  with  (a)  arbitrary 
non-stationary  Gaussian  random  signals,  (b)  Gaussian  Markov  meshes,  (c)  non-Gaussian 
Markov  meshes. 

All  these  results  were  obtained  under  the  hypothesis  of  the  complete  knowledge  of 
the  probability  distnbutions  of  the  random  signal  being  approximated.  In  the  case  when 
the  random  process  is  given  with  the  help  of  correlations  and  cross-correlations,  we  face 
a  new  problem  the  solution  of  which  requires  different  methods.  At  present,  it  seems  that 
the  general  case  of  arbitrary  processes  of  this  kind  is  intractable,  except  in  some  particular 
cases  such  as  the  Gaussian  processes.  In  those  cases  we  are  particularly  interested  in 
harmonizable  processes  that  have  well-defined  joint  time-frequency  representations,  and 
this  research  would  investigate  the  new  approximation  problem  for  such  processes. 

During  the  last  year,  specifically,  M.  Rosenblatt-Roth  obtained  a  senes  of  important 
results  in  these  areas  as  well  as  in  some  new  directions  of  research.  In  particular, 
concerning:  (a)  the  concept  of  relative  entropy  convergence;  (b)  the  concept  of  a 
continuous  projector;  (c)  the  best  approximation  by  means  of  Gaussian  probability 
measures;  (d)  the  study  of  the  structure  of  the  Boolean  algebra  of  Markov  mesh  projectors; 
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(e)  the  study  of  the  structure  of  the  Boolean  algebra  of  Gauss-Markov  mesh  projectors;  (f) 
the  amount  of  Gauss-Markov  mesh  information;  (gj  _  ^pt  of  decreasing  stochastic 
dependence;  (h)  the  concept  of  intensity  of  stochastic  dependence;  (1)  clustering  of  random 
variables;  (j)  conditional  information;  (k)  the  general  explicit  expression  of  a  Gauss- 
Markov  mesh. 

In  addition,  M.  Rosenblatt- Roth  defined  and  began  to  study  of;  (a)  the  important  class 
of  completely  regular  projectors;  (b)  the  best  approximation  of  probability  measures  with 
various  classes  of  Markov  meshes;  (c)  the  extension  of  the  theory  of  Markov  meshes  to 
the  case  when  when  the  characteristic  sets  and  kernels  are  random;  (d)  the  extension  of 
the  best  approximation  method  to  generalized  random  processes  and  to  random  processes 
with  continuous  time;  (e)  the  best  approximation  of  noisy  channels. 
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